{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\r\n",
    "from numpy.linalg import det\r\n",
    "\r\n",
    "from functools import reduce\r\n",
    "import MatDisp "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def LU(A:np.ndarray):\n",
    "    assert A.ndim == 2\n",
    "    n = A.shape[0]\n",
    "    assert n == A.shape[1]\n",
    "    L = np.eye(n,n)\n",
    "    U = np.zeros((n,n))\n",
    "    for r in range(n):\n",
    "        U[r,r:n] = A[r,r:n] - L[r,0:r]@U[0:r,r:n]\n",
    "        L[r+1:n,r] = (A[r+1:n,r] -L[r+1:n,0:r]@U[0:r,r])/U[r,r]\n",
    "    return L,U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": "$$\\begin{bmatrix}10&-7&0&1\\\\-3&2&6&2\\\\5&-1&5&-1\\\\2&1&0&2\\end{bmatrix}$$",
      "text/plain": "array([[10.      , -7.      ,  0.      ,  1.      ],\n       [-3.      ,  2.099999,  6.      ,  2.      ],\n       [ 5.      , -1.      ,  5.      , -1.      ],\n       [ 2.      ,  1.      ,  0.      ,  2.      ]])"
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([\n",
    "        [10,-7,0,1],\n",
    "        [-3,2.099999,6,2],\n",
    "        [5,-1,5,-1],\n",
    "        [2,1,0,2]\n",
    "    ])\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "L,U=LU(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": "$$\\begin{bmatrix}1&0&0&0\\\\0&1&0&0\\\\0&-2500000&1&0\\\\0&-2400000&1&1\\end{bmatrix}$$",
      "text/plain": "array([[ 1.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00],\n       [-3.0000000e-01,  1.0000000e+00,  0.0000000e+00,  0.0000000e+00],\n       [ 5.0000000e-01, -2.5000000e+06,  1.0000000e+00,  0.0000000e+00],\n       [ 2.0000000e-01, -2.4000000e+06,  9.5999968e-01,  1.0000000e+00]])"
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": "$$\\begin{bmatrix}10&-7&0&1\\\\0&0&6&2\\\\0&0&15000005&5749998\\\\0&0&0&5\\end{bmatrix}$$",
      "text/plain": "array([[ 1.00000000e+01, -7.00000000e+00,  0.00000000e+00,\n         1.00000000e+00],\n       [ 0.00000000e+00, -1.00000000e-06,  6.00000000e+00,\n         2.30000000e+00],\n       [ 0.00000000e+00,  0.00000000e+00,  1.50000050e+07,\n         5.74999850e+06],\n       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n         5.07999891e+00]])"
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U"
   ]
  }
 ],
 "metadata": {
  "language_info": {},
  "orig_nbformat": 3
 },
 "nbformat": 4,
 "nbformat_minor": 5
}